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1. INTRODUCTION 

In power systems, both active and reactive power demands are never steady and continuously 
change with rising or falling trend. Steam input to turbo generators (or water input to hydro generators) must 
therefore, be continuously regulated to match the active power demand, failing which the machine speed will 
vary with consequent change in frequency, which may be highly undesirable. In brief, the changes in real 
power affect the system frequency, while reactive power is less sensitive to changes in frequency and is 
mainly dependent on changes in voltage magnitude. The quality of power supply must meet certain minimum 
standards with regard to constancy of voltage and frequency [1]. 

The operational objective of LFC is to maintain reasonably uniform frequency, to divide load 
between generators and to control the tie-line interchange schedules. The change in the frequency and tie-line 
power are sensed, which is a measure of the change in rotor angle, i.e., the error to be corrected. The error 
signal, i.e., Af and APtie, are amplified, mixed, and transformed into a real power command signal APv, 
which is sent to the prime mover to call for an increment in the torque. The prime mover, therefore, brings 
change in the generator output by an amount APg which will change the values of Af and APtie within the 
specified tolerance. 

In this paper mathematical model for the two area system is developed in order to analyses and 
design the control system. A state space model is developed by linearizing the mathematical equations of 
different components and thus the state matrix is formed. An Eigen value based objective function to place 
the modes in better region on the complex plane is considered and optimum controller parameters are found. 
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2. SYSTEM INVESTIGATED 
During normal operation, the real power transferred over the tie-line is given by (as shown in Figure 
1 and 2), 


|Z, ||E.| . 
P a= sind (1) 
12 X: 12 





Where, ô =ð -ó 
For a small deviation in the tie-line flow from the nominal value: 


dP 
AP, -( #2 Jaa, (2) 
= PA, 





Where Ps is the slope of power angle curve at initial operating angle called synchronizing power coefficient. 


E,||E,| 





Fy = | X, cos AG, (3) 
AF, = P,(6,—6,) (4) 
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Figure 2. Electrical equivalent for a two area system 


A block diagram representation of two area system with LFC containing only primary loop is shown 
in Figure 3 with each area represented by an equivalent inertia M, load damping constant D, turbine and 
governing system with an effective speed droop R. Tie-line is represented by synchronizing torque 
coefficient. A positive represents an increase in power transfer from area 1 to area 2 and it is equivalent to 
increasing load in area 1 and decreasing load in area 2. 
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Figure 3. Two area system with primary LFC loop 


Both areas will have a same steady state frequency deviation for a load change of in area 1. 





Aw = Ao — 4o, and AP — AP, —AP,, =A@D, (5) 


AP,, + AP, = A@D, (6) 





From the governor speed characteristics the change in mechanical power is: 


AP,, =Aa/ R, (7) 
AP,,, =Aa/ R, 

8 
hee AP, (8) 








Where, 
1 
a-(f+0) (9) 
1 
B, -(4+0,] 
l D, | AP, 
R, 2 LI (10) 
AP 








In normal operating condition power system is operated so that demand of areas is satisfied at 
nominal frequency. A simple control strategy should include the following functions: 
- Frequency approximately at nominal value. 
- Maintaining the tie-line flow at about schedule. 
- Each area should absorb its own load changes. 


2.1. Tie-line bias control 

Conventional LFCs are based on tie-line bias control where each area tends to control its area 
control error (ACE) to zero. The control error consists of linear combination of frequency and tie-line 
error [2]. 
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ACE, = >"._, AP, + K,A@ (11) 


The area bias factor K; determines the amount of interaction during the disturbances in the 
neighboring area. In order to get satisfactory performance area bias factor is selected as: 


RHC ep) (12) 
R, 

ACE, = AP, + BAO, l 

ACE, = AP, + B,Aa, (13) 

AP, and AP,, are the deviations from the scheduled interchanges. ACEs are the actuating signals that are 


used to change the reference set points, and when the steady state is reached and will be zero. The block 
diagram of two area system LFC with integral control is shown in Figure 4. 
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Figure 4. Block diagram of two area system LFC with integral control action 


In this paper the integral control block is replaced by PID controller and state space model of the 
system X=AX is obtained. State matrix of order 11x11 is formed with state variables AP AP, AP 


ml?’ 





Aa > ACE,> AP, AP, AP, AP p > A,» ACE, respectively. The parameters of the system which are 
to be optimized using PSO are K, Kno Kp Re By Kp> Ko» Kpo> Ry? B, respectively [3]. Block 


diagram of two area system with PID controller is shown in Figure 5. 
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Figure 5. Block diagram of two area system LFC with PID control action 





K 
PID1(s) =K,, +— + Kps (14) 
Ss 
K 
K a ++ Kps 
PID2(s) = S (15) 


3. PARTICLE SWARM OPTIMIZATION 

The PSO algorithm was first introduced by Dr. Russel C. Eberhart and Dr. James Kennedy (1995), 
inspired by social behavior of bird flocking or fish schooling. The system is initialized with a population of 
random solutions and searches for optima by updating generations [4]. However, unlike GA, PSO has no 
evolution operators such as crossover and mutation. In PSO, the potential solutions, called particles, fly 
through the problem space by following the current optimum particles. Each particle keeps track of its 
coordinates in the problem space which are associated with the best solution (fitness) it has achieved so far. 
(The fitness value is also stored) this value is called pbest. Another "best" value that is tracked by the particle 
swarm optimizer is the best value, obtained so far by any particle in the neighbors of the particle. This 
location is called Ibest. When a particle takes all the population as its topological neighbors, the best value is 
a global best and is called gbest. 

In the numerical implementation of this simplified social model, each particle has three attributes; 
the position vector in the search space, the current direction vector, the best position in its track and the best 
position of the swarm. The process can be outlined as follows [5]: 

Step 1. Generate the initial swarm involving N particles at random. 

Step 2. Calculate the new direction vector for each particle based on its attributes. 

Step 3. Calculate the new search position of each particle from the current search position and its new 
direction vector. 

Step 4. If termination condition is satisfied, stop. Otherwise, go to step 2. 

As the particle can fly in D-dimension search space, the position and velocity of i-th particle can be 
represented as: 


Xi=[Xi1,Xi2,Xi3,Xia,....... xiD] (16) 
Vi=[vil,vi2,vi3,vi4,...... viD] (17) 


With increased iteration, the swarm will move towards its global best position by keeping track of 
their personal best. In D- dimensional search space the pbest of i-th particle can be represented as: 


Phest=[Dpit,Pi2,Pi3,Pi4, - aie . pip] (18) 
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and gbest of the whole swarm is presented as: 


gbest=[81,82,83; Z4... 0.000 gp] (19) 


The new direction vector of the i-th particle at time t, is calculated by the following scheme 
introduced by Shi and Eberhart. 


Vi" = œv, +clRi(pbesti, — x4) eo 
+c2R} (gbest, — x) 


Ri and R; are random numbers between 0 and 1. vi and PA is the velocity and position of the i-th 
particle in d-th dimension at its time track t. pbesr:, 18 the best position of the i-th particle (personal best) in 
d-th dimension in its track at time t and gbest', is the best position of the swarm in d-th dimension at time t. 


There are three parameters such as the inertia of the particle œ , and two parameters cl and c2. c1 and c2 are 
the learning factors which determines the relative influence of the cognitive and social components to update 
the position and velocity component. 

Then, new position of the i-th particle at time t, Ba Scie is calculated from: 


Xa =X a eV, 21) 


Where X/, is the current position of the i-th particle at time t. After the i-th particle calculates the next 
search direction vector yo in consideration of the current search direction vector vi, othe direction vector 
going from the current search position X; to the best search position in its track pbest', and the direction 
vector going from the current search position Xj, to the best search position of the swarm gbest', » it 
moves from the current position x+, to the next search position X/}! calculated by Equation (21). In 


general the parameter œ is set to large values in the early stage for global search, while it is set to a small 
value in the last stage for local search. 

The inertia weight is used to control the impact of the previous velocities on the current velocity, 
influencing the trade-off between the global and local experience. Although Zheng claimed that PSO with 
increasing inertia weight performs better, linear decreasing of the inertia weight is recommended by Shi and 
Eberhart. 


Iter vax 


Fae [a= ) Pee (22) 


Where W wax and W „in Ve Maximum and minimum of inertia weight value respectively, iter,,,, 18 Maximum 


iteration number and iter is the current iteration. A so-called constriction factor K, is factor that increases the 
algorithm’s ability to converge to a good solution and can generate higher quality solution than the 
conventional PSO approach. In this case, the expression used to update the particle’s velocity becomes 





Vij! = K *[o'vi, + c1R! (pbest!, — x4) (23) 
+c2R)(gbest', — x}, )] 
Where 
K= 2 I 
2-9-(¢" ap)! 
p=cl+c2, g>4 = 
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4. PROBLEM FORMULATION 

Provide a statement that what is expected, as stated in the "Introduction" chapter can ultimately 
result in "Results and Discussion" chapter, so there is compatibility. Moreover, it can also be added the 
prospect of the development of research results and application prospects of further studies into the next 
(based on result and discussion). 


4.1. Objective function 
According to appendix B the equation of D-contour is given by 


F(z) =Re(z)—min[-C Img(z) |, a] = 0 (25) 


Where z € A is a point on the D-contour and A represents the complex plane. 
Defining J as: 


J=max[Re( 7 )-min(-¢lIm 4 |,a)] (26) 


Where n is the number of Eigen values. 4, is the i-th Eigen value of the system at an operating point. A 


negative value of J implies that all the Eigen values lie on the left of the D-contour. If J is positive that 
implies Eigen value is lying on the right of contour. 
On these facts objective function F can be defined as: 


F=-|J if J<O (27) 
a if J>O 


Where a is a large positive number. 
The optimization problem can now be stated as: 
Minimize F Subject to, 


min max 

K <K„<K” 
min max 

Ki SK, SK; (28) 
min max 

K pi < Kp, = Kp; 


Re < R, < Ree 
B, =— +D, (29) 


Where Kyi, Ki, Kpi are PID controller parameters Ri and Bi are speed characteristics and area bias factor 
respectively. i=1, 2 for control area first and second respectively. 


4.2. Algorithm 

Stepl. Initialize the set of particles with position value (population), its upper limit and lower limit, random 
velocities, inertia weight, acceleration constants, itermax, etc. 

Step2. Linearize system and calculate the eigen values for each particle from the system modelling. 

Step3. Calculate the objective function J. 

Step4. Set all position values as local best values and fitness values of objective function as local fitness. 
Find global fitness and its corresponding position value. 

Step5. Set iter=1 and compute inertia weight by 


ae (1a ) ees (30) 


Iter vax 


Update velocity using velocity update equation for all particles by 
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Step6. 


Step7. 


Step8. 


Step9. 


Vi" = o'vi, +clRi (pbesti, — x) G1) 
+c2R}(gbest', —xi,) 

Also check its upper and lower limits. 

Update position value of particles by 

Xa = Xia + Via" (32) 


Also check its upper and lower limits. 

Linearize system and calculate Eigen values for each particle and corresponding fitness values of 
objective function. 

Now compare these fitness values to the previous fitness values. Minimum fitness values will be 
selected as local fitness values and its corresponding position values as local best values. Find 
minimum of all fitnesses, e.g. global fitness value and its corresponding position value. 

If number of iteration reaches itermax go to step 10, otherwise go to step 5. 


Step10. Particle with minimum fitness value is the optimum particle. 


4.3. Result 


PID parameters and speed characteristics of both areas are allowed to vary within following ranges: 


-5.00 < K, < 5.00 

—5.00 < K, < 5.00 (33) 
—5.00 < K, <5.00 

0.02 < R < 0.08 


Time constants of different model components associated with both areas are given in appendix A. 


The desired PSO parameters are given in Table 1. with the desired setting of damping ratio ‘6’ (0.5) 


and reference line value ‘a’ (-1.5) algorithm is run several times. Among all the runs most optimum run and 
its corresponding parameters are selected. Optimized parameters and corresponding Eigen values as shown in 
Table 1-3. 


Table 1. Optimized parameters 





Ka -1.8787 
Ky -2.3640 
Kopi -0.9364 
Ri 0.0652 
Bı 15.9301 
K -0.1353 
Kp -0.5873 
Kaz -0.5212 
R2 0.0653 
B2 16.2087 





Table 2. Optimized Eigen values (oscillatory mode) 
Eigenvalues Damping 
-1.4953 + 2.41921, -1.4953 - 2.41921 0.5258 
-1.4934 + 1.35701, -1.4934 - 1.3570i 0.7401 
-1.4762 + 0.85551, -1.4762 - 0.8555i 0.8652 
-1.4887 + 0.5015i, -1.4887 - 0.5015i 0.9477 
-0.2653, -0.0600, -0.1125 1 
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Table 3. PSO parmeters 
CI 2 
C2 2 
Wmax 0.90 
Wmin 0.10 
Vmax 100% of Pmax 
Vmin 100% of pmin 
Population 500 
Iterations 1000 
CI 2 





5. SIMULATION STUDY 
The system is simulated in MATLAB environment for step change in load (0.2 pu) in area first, 


given at one second and is checked out for 50 seconds. Comparison of responses of Aw,, AP, ACEL, AP 


,Aw,, ACE2, AP, between conventional integral control and PIDPSO is carried out as shown in 


Figure 6-12. 




















Figure 6. Frequency deviation response of area first with integral control (dashed line) and PID control 
(solid line) 














Figure 7. Power interchange response of area first with integral control (dashed line) and PID control 
(solid line) 
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Figure 8. Area control error (ACE) response of area first with integral control (dashed line) and PID control 
(solid line) 





Figure 9. Deviation in mechanical power response of area first with integral control (dashed line) and PID 
control (solid line) 





Figure 10. Frequency deviation response of area second with integral control (dashed line) and PID control 
(solid line) 
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Figure 11. Area control error (ACE) response of area second with integral control (dashed line) and PID 
control (solid line) 





Figure 12. Deviation in mechanical power response of area second with integral control (dashed line) and 
PID control (solid line) 


6. CONCLUSION AND FUTURE WORK 

The target of the developed work is to improve the control performance of two area power system 
using a controller based on evolutionary algorithm. PIDPSO shows the better control performance than 
conventional integral control in terms of settling time, over/under shoots. The discussed formulation can be 
extended to multiarea power systems. 


APPENDIX A 
From the block diagram shown in Figure 5. the state space equations can be written as: 
1 
AP, (s)=————- AP (s (34) 
x (s) 1+T,(s) 6s) 
AP, (s)=—-— AP, (s) 
ig 1+ T, Cs) a (35) 
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1 
(AP, (s) -AR (s) —— _— = Als) 
i (2H +D) (36) 
Aa(s)B, —AP,,(s) = ACE(s) (37) 
K, 1 
ACE(s)(K, ++ K,)-— A@(s) = AP, (s) 
s R : (38) 
System data- as shown in Table 4. 
Table 4. System parameters 
Tgi 0.2 sec 
Tez 0.3 sec 
Tu 0.5 sec 
Te 0.6 sec 
H1 5 
H2 4 
D1 0.6 
D2 0.9 
Base power 1000 MVA 
Ps (synchronizing power 2.0 pu 
coefficient) 
APPENDIX B 
From Figure 13. the equations of line ABM and DCM in x-y plane can be written as: 
aa (39) 
x=- 
a: | y| +x=0 
k (40) 
Where is the modulus of the slope of line. The equation of the line BC can be written as: 
x-a=0 (41) 
Combining Equations (39) and (40) the equation of D-contour ABCD can be written as: 
x—min(- ly I/k, a)=0 (42) 
If -1/k=C, then in reference of complex plane the equation (41) becomes: 
Re(z)—min(-C IImg(z) |, a)=0 (43) 


Where x=Re(z), and y=Img(z). z=x+iy 
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x -axis 
Figure 13. D-contour in x-y plane 
Dii;63. 
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